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Abstract 



We consider two strongly hyperbolic Hamiltonian formulations of gen- 
eral relativity and their numerical integration with a free and a partially 
5-H , constrained symplectic integrator. In those formulations we use hyper- 

bolic drivers for the shift and in one case also for the densitized lapse. 
A system where the densitized lapse is an external field allows to enforce 
■ the momentum constraints in a holonomically constrained Hamiltonian 

^ ' system and to turn the Hamilton constraint function from a weak to a 

0^ , strong invariant. 

. These schemes are tested in a perturbed Minkowski and the Schwarz- 

schild space-time. In those examples we find advantages of the strongly 
hyperbolic formulations over the ADM system presented in [28]. Further- 
^SJ ' more we observe stabilizing effects of the partially constrained evolution 

' in Schwarzschild space-time as long as the momentum constraints are en- 

. forced. 



PACS numbers: 04.25.D-, 04.20.Fy 



■ 1 Introduction 



Many of the recently developed schemes in numerical relativity are based on the 
BSSN system 29, ,5^. One of the key features that distinguishes BSSN from e.g. 
the dynamical ADM equations is its strong hyperbolicity. Systems with that 
property are preferred, because the wcU-posedness of the Cauchy problem can 
be proven [25j and often they lead to numerical schemes with improved stability 
properties [T3j . Using strongly hyperbolic systems stable evolutions of the black 
hole binary problem are possible since 2005 pSl [Til lil [T^ . 

Here we are interested in numerical schemes that apply symplectic integrators 
for Hamiltonian systems. These integration methods respect the Hamiltonian 
structure in the discretization and have been essential for attaining favorable 
propagation properties in various areas of scientific computing (see, e.g., [211 
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[531 HZ] and references). The problem there is that, m contrast to the dynamical 
ADM equations |31 which are just weakly hyperbolic, many strongly hyperbolic 
formulations known do not possess a Hamiltonian structure. Yet, there is one 
exception, namely the Hamiltonian generalized harmonic system that Brown 
proposed in [TT] . 

In this article we present a second strongly hyperbolic Hamiltonian system, 
the fixed lapse system, and apply the free and constrained symplectic integra- 
tors introduced in to these two systems. Other applications of symplectic 
integrators in numerical relativity are presented e.g. in [71 [51 ITOl [T7] . 

We compare the properties of our schemes in numerical experiments with 
effectively l+l-dimensional versions of Einstein's equations. The purposes of 
these experiments are the following. At first we want to confirm the expecta- 
tion that schemes based on the strongly hyperbolic systems lead to more stable 
results than schemes based on the (just weakly hyperbolic) ADM equations. 
Furthermore we want to recover the good propagation properties of nearly con- 
served quantities that we found in [28j with the new schemes here. Finally, the 
major question addressed is, for which systems and in which situations con- 
strained evolution is beneficial in comparison with free evolution. In |28j we 
found better results for constrained integration of the ADM equations, but it is 
not clear whether this is still valid when the free evolution schemes are based 
on strongly hyperbolic systems. 

We start our considerations in section 2 where we introduce the two Hamilto- 
nian formulations and discuss their pros and cons. Brown's Hamiltonian has the 
advantage that it deals with well-known gauge conditions, whereas for the fixed 
lapse system we can apply a numerical scheme that enforces the momentum 
constraints in a holonomically constrained Hamiltonian system. 

In section 3 we describe the spatially discrete Hamiltonians and the sym- 
plectic integration methods used. In the free evolution schemes the integrator is 
the Stdrmer-Verlet method, a standard second order symplectic integrator (see 
[501 HH 121] ) • There the Hamiltonian and momentum constraints are not taken 
into account and may drift off. 

The problem of constraint growth is addressed by applying the RATTLE 
method, a constrained version of the Stormer-Verlet integrator. It is applicable 
to holonomically constrained Hamiltonian systems. For Brown's Hamiltonian 
it can be used to enforce a combination of the momentum constraints and a 
set of new constraints that arise since the densitized lapse becomes a dynamical 
variable. For the fixed lapse system we can achieve even more, namely that the 
momentum constraints themselves are enforced. The mechanism behind is the 
same as for the ADM equations [25]. 

In section 4 we discuss the important aspects of the l+l-dimensional imple- 
mentation in detail and introduce the test examples with which we have done 
our numerical experiments: a perturbed Minkowski problem and Schwarzschild 
space-time in two different coordinate systems. 

Section 5 compares four integration methods for a perturbed Minkowski 
problem. It turns out that in three of the four schemes one obtains similar 
behavior, only constrained evolution of the fixed lapse system leads to instabili- 
ties that presumably can be explained through an inappropriate choice of initial 
data. 

In sections 6 and 7 we compare the free and the constrained symplectic 
schemes on the Schwarzschild space-time. Section 6 deals with Brown's gener- 
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alized harmonic system, showing that the resuhs of both schemes are similar. 
In section 7 the fixed lapse system is considered. There it turns out that er- 
rors which arise at the boundaries in the free evolution are suppressed in the 
constrained scheme. 

Our numerical experiments thus confirm this remarkable property of the 
constrained scheme that enforces the momentum constraints. But they also 
show that constrained evolution is not beneficial when applied to a constrained 
Hamiltonian system that does not enforce the momentum constraints. 

We finish the article by showing how the fixed lapse Hamiltonian can be 
derived in appendix A. 

2 Hamiltonian formulations of general relativity 
with dynamical shift and densitized lapse 

As motivated in the introduction we are interested in Hamiltonian formulations 
of general relativity that provide strongly hyperbolic equations of motion. There 
is no unique way to find such a formulation. We follow the ideas proposed by 
Brown [TT]. That is, we include hyperbolic drivers for the shift and the densi- 
tized lapse functions. Our approach differs in some details, and we summarize 
it in what follows. 

The systems of interest are second order in space and first order in time. 
A definition of strong hyperbolicity for those systems is given e.g. in [19j . For 
strongly hyperbolic systems it can be proven that they possess a well-posed 
Cauchy problem |251 113j. they are hence preferred in numerical applications. 

In the construction we start from the well-known Super- Hamiltonian \15\ 116] 
(see also [21 HI])- There the position variables are provided by the 3-metric hij 
and the corresponding canonical momenta are denoted tt*^ . They are related 
to the extrinsic curvature Kij through tt*-' — \fh {K^iW^ — if*-'). Here h is the 
determinant of the 3-metric hij. The Super-Hamiltonian takes the form 

T^s (aC + /3'a)rf^a; (1) 

with freely specifiable (external) functions a and The densitized lapse a is 
related to the lapse function N through a = N/\/h, and is the shift vector. 
The functions C and Ci are given by 

C = TT''^^,-]^lT\n^,-hR, (2) 

where R is the Ricci scalar of the metric hij , anc(3 

Q = -2h,jDkTT'^. (3) 

The equations C = and C,; = are the Hamiltonian and momentum con- 
straints respectively. 

-"^ Notice that tt-'* is a tensor density of weight +1. 
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2.1 Shift and densitized lapse as dynamical variables 

According to [TT] we can extend the phase space and use a modified Hamilto- 
nian such that shift and densitized lapse become dynamical variables. Here we 
treat the shift as a momentum and the densitized lapse as a positionH We intro- 
duce their canonically conjugate variables 7^ and a respectively. The canonical 
symplectic two-form for the extended phase-space of (hij, a, 7^; tt*-' , ct, is then 
dhij A d-Tr'^ + da f\ da + d-ji A d/3^ . 

The interpretation of /3* as momentum variables may lead to complications 
if one wants to transform from the Hamiltonian equations of motion to an 
action with Euler-Lagrange equations. The 3-1-1 decomposition of the Einstein- 
Hilbert action explicitly depends on the shift. As an action it must not contain 
momentum variables, but only positions and their time derivatives. Yet, given a 
Hamiltonian it depends on its concrete form whether the shift can be expressed 
in terms of positions and their time derivatives. Hence, it may be impossible to 
perform the Legendre transformation to obtain an appropriate action. 

Here we take the Hamiltonian formulation of the equations as fundamental 
and do not consider possible problems concerning the action and corresponding 
Euler-Lagrange equations further. 

Yet, in [TT] Brown gets the constraints 

cr = and 7, = (4) 

as primary constraints from the Einstein-Hilbert action. Since for us the Hamil- 
tonian formulation is fundamental we are not limited to this case. However, of 
course we are interested in an Hamiltonian that provides dynamical equations 
for hij and tt*^ that are consistent with the equations of motion in general rela- 
tivity. We will see that if ^ is satisfied then this requirement is automatically 
fulfilled with an appropriate extended Hamiltonian. 
We define 

Hg := J d^x (Ad + f]*7,) , (5) 

where A, fi* are functions of the canonical variables that are linear in tt*^ , cr, 7^ 
and in the first spatial derivatives dkhij, dia, di(3^ ^ with coefficients that depend 
on hij, a and /?■'. 

Due to the structure of the gauge Hamiltonian Tig the canonical Hamiltonian 
equations of motion for hij and tt*-' that correspond to the full Hamiltonian 
Tis + Ti-g are equivalent to the dynamical part of the ADM equations when the 
constraints ([U are satisfied. But the freedom in Hg can be used to introduce 
quite general hyperbolic drivers for a and 

For concrete calculations we need to make one particular choice for the func- 
tions A and . The simulations that we discuss here are based on a generalized 
harmonic system. We set 

A = - ad^P' + ^a^a (6) 

and 

= ~P^dj(3' - ha^T)kh^^ + hah'Wja - ^ha^h'^-f^. (7) 

^The reasons for this choice are discussed in section [2. 1.1 1 or in I28| . 
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There are also other possible choices, and in Brown discusses some of them. 

The canonical Hamiltonian equations of motion, i.e. the generalized har- 
monic system, is now 

hij = 2aTTtj - an'^khij + Di/3j + DjP^, (8a) 
+ hD'D^a - hh'W^Dua + hDiiP^n'^) 

1 . 1 n 1 1 1 

— -a h^J^ + -a hh'-^^k^ — -ahj^d^a — -ahYd^a 

+ a^hh'^h'^'dkli - ^a^hh''h^"'j''dihkm - ^a^hh''h^"'-f''dmhki, (8b) 
a = a^a + P'^,a-a^^P\ (8c) 
& = -C + P'd.a + 2(jd,l3' + ahh'^d^-yi 

- ^ (aV^ - a^h-tajh'^) - ahy^h'^-fid^hjk + ah-fih'^ h^'' d^hjk, (8d) 
7, = + 2(Td,a + ad,<j - 7,a,/3^' + j,djp' + p^d.j,, (8e) 
j3' = l3^dj/3' - a^hVi^^h'^ + a^hT]t,h^^ + a^hh'^jj - ahh'Wja. (8f) 



We see that ([8a|) is also one of the ADM equations, i.e. the dynamical behavior 
of hij is not changed. The equation (jSbp for tt*-' differs from the corresponding 
ADM equation, it involves additional terms. But these terms vanish when 7^ 
vanishes. Hence, in order to describe the dynamical behavior of tt*-' correctly 
we need to require that 7^ = 0. 

In (ID) we not only have the constraint 7^ = 0, but also u = 0. Yet, the 
equations of motion for the "physical" variables hij and tt*^ already agree with 
the ADM equations when only 7,; = is satisfied. One may hence conclude that 
a does not need to vanish. 

However, for cr 7^ we find instabilities already when hij and tt'-' are exact 
solutions of the ADM equations. In particular, when 7^ = and the vector 
constraints are satisfied then equation (j8ep implies 



= 2ad^a + ad^a. (9) 

It follows that (t(x, i) = c{t) / {x, t) . Then, if the scalar constraint is satisfied 
as well we obtain from (I5cl) and (I8dl) that 



'^(*) = o "77 — TT' where cq = c{to). (10) 

2 - co(t - to) 

Hence, we obtain stable results only when |co(i — io)| remains small. Therefore 
it is obvious that ct = should be satisfied, too. The full set of constraints to 
be satisfied is thus 

C = 0, a = 0, (7 = 0, 7» = 0. (11) 
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2.1.1 Constrained integration of the generalized harmonic system 

In [2S] we proposed a numerical integration scheme for the ADM equations that 
enforces the momentum constraints in a holonomicaUy constrained Hamilto- 
nian system and turns the Hamilton constraint function from a weak to strong 
invariant. 

The main steps to achieve this were to interpret the shift as a momentum 
and to impose the holonomic constraints 7^ = as well as an appropriate gauge 
condition. The obvious question is now whether an analogous scheme with 
similar properties can be constructed also for the generalized harmonic system 
©• 

Unfortunately this is not the case. For the ADM equations the momentum 
constraints could be enforced, because they were hidden constraints of 7^ = 
(i.e. 7i = Ci). Here we do not have this direct relation, but we have ([5e|l instead. 
Hence, when we enforce 7^ = at each time then we only get 

= a + 2(T^^a + ad^a, (12) 

i.e. only this combination of the momentum constraints and cr = can be 
enforced. This might still be beneficial, but we cannot control the momen- 
tum constraints separately and we cannot expect that the Hamilton constraint 
function becomes a strong invariant. 

Obviously the generalized harmonic system ((HJ here behaves differently than 
the ADM equations because of the terms that contain a and cr. These terms will 
virtually always appear when the densitized lapse is treated as a dynamical vari- 
able. Therefore we investigate in the next section a Hamiltonian system where 
a is an external field, and only /3* is subject to a dynamical gauge condition. 



2.2 Hyperbolic shift driver with fixed densitized lapse 

Since we cannot enforce the momentum constraints in a holonomicaUy con- 
strained Hamiltonian system when the densitized lapse function is a dynamical 
variable, we are interested in a Hamiltonian formulation where a is an exter- 
nal field. Thus, we consider the phase space of (hij,"fi]Tr^^ , (3^) with syniplectic 
two- form dhij A dn^^ + d'ji A . 

We define another gauge Hamiltonian 

= y d^xd'-f,, (13) 

and need to find a function such that the canonical Hamiltonian equations 
of motion that correspond to Tis + "H/s become strongly hyperbolic. 

It turns out that this requirement is very restrictive when a is an external 
field. We are able to construct an appropriate function: 

= -p^d.P^ - ^a^T)^hh=^ + ^-a^T^hh'^ (14) 

2 2 
-I- -ahh'Wja + -a^hfi'^jj, 

but every other possible choice for that we investigated leads to a system 
with the same principal part as (jlSp (see appendix IX]) . 
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The canonical Hamiltonian equations of motion that we derive from TLs+Hp 
using are now 

hij = 2a'rTij - an'^khij + Di(3j + DjPi, (15a) 
7T*^' = 2a7r'fc7r^'= - air''^ t:^ k + ahB!^ - a/i/i'^i? (15b) 
+ T^^^DkP' + tt'^DuP^ - Tr'WklS'' - P'^Dk-n-'^ ~ hD'D^a + hh^^D^Dua 

+ y (^ahTi^{h^^h^'^ - /i'^7i'™)7fc - 3a/ir|„(/i""7^' + /i^"7* - 2/1*^7") 

- 2^2/17*7^ + 2a^hh'^^^-ik + 2ahD'-f^ + 2ahD^f - Gahh'Wkj'' 

7» = a- 7,9./5^' + 7»5j/3'' + f3'd,j,, (15c) 
= -/3J9j/3^ + ^ (Sa^r^:^^/!/!^'^ - 2a^T)j,hy^ + 2a^hj' + ahD'a^ . (15d) 

In what follows we denote these equations the fixed lapse system. Again (|15ap is 
one of the ADM equations and (jl5b[) becomes an ADM equation when 7^ = 0. 

We also see from ()15cp that the momentum constraints vanish when 7^ van- 
ishes identically. Therefore, in contrast to the generalized harmonic system (jS]) 
we can enforce the momentum constraints in a holonomically constrained sys- 
tem. Since the dynamical behavior of hij and tt'^ is then equivalent to the ADM 
equations, we also get that the Hamilton constraint function becomes a strong 
invariant [2]: 

[dt - l3'Di)C = 0. (16) 

Even if this property does not extend to the space discretization, it is an extra 
bonus for this momentum-constrained formulation. 

In the numerical examples we will use 1+1 dimensional simplifications of the 
two formulations ^ and (fT5|l . 

3 Discrete Hamiltonian and numerical integra- 
tion methods 

To apply numerical methods one approximates the continuous functions hjj^ 7^, 
a, TT*-' , and a through objects with finitely many degrees of freedom^ We 
collect the finite number of unknowns in four vectors q, p, 7 and (3. For the 
discretized system (|S]) q and p correspond to the discretizations of (hij , a) and 
(tt*-' , a) respectively, whereas for the system (IT5t these vectors correspond to 
the discretizations of hij and tt*-' respectively. The ordering in q, p, 7 and /3 is 
chosen such that components corresponding to the same grid point are ordered 
consecutively. 

For both considered systems the full Hamiltonian {Jis + Tig respectively 
Tis+Ti-p) consists of terms that are either quadratic in the momentum variables 
(tt*-' , cr, /3') or independent of them. Any reasonable discretization of the full 

^We use finite difl^erences and piecewise constant functions, as described in section [4. II 
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Hamiltonian will hence assume the forir0 

^(q, P) = ^P^S(q)p + C/(q, 7) + /3^D(q)p + i/3^E(q, j)f3, (17) 

where S(q), D(q) and E(q, 7) are matrices of the appropriate dimensions. 
The canonical equations of motion for this discrete Hamiltonian are 

q=S(q)p + D(q)^/3, (18a) 
7 = E(q,7)/3 + D(q)p, (18b) 

P = -^P^VqS(q)p - VqC/(q,7) - /3^VqD(q)p - i/3^VqE(q,7)/3, (IBc) 



/3=--/3^V^E(q,7)/3- V^C/(q,7). (18d) 



In section 13.11 we introduce a free symplectic method that integrates ^ 
numerically and in section 13.21 a constrained scheme is presented. The latter 
additionally requires the following constraints 

7 = 0, g(q) = 0. (19) 

The equation g(q) = is a gauge condition and fixes /3. It must be chosen such 
that the following matrix is invertible (see also [HI [H] ) 

A(q) := Vqg(q)^D(q)^. (20) 

Time differentiation of g(q) = and using p8ap for q gives 

Vqg(q)^S(q)p + A(q)/3 = 0, (21) 

which shows that indeed (3 is determined by the gauge condition, when A(q) is 
invertible. 

A candidate for the choice of the function g is a discretization of the Dirac 
gauge, dj{h^/^h^i) = [Ml [SI- With this choice, A(q) is a discretized second- 
order elliptic differential operator. In section [4. II we discuss other gauge condi- 
tions with a similar structure for a simplified system. 

3.1 The Stormer-Verlet method 

A standard symplectic integrator for Hamiltonian systems is the Stormer-Verlet 
scheme (see, e.g., [20 ). When applied to (fTS]) . a step from values (q", 7", p", /3") 
at time T to values (q"+\ 7"+\ p"+\ /3"+^) at time = + At reads as 
follows: 

p„+i/2 ^ p"-^(i(p"+i/2)^VqS(q")p"+i/2 + VqC/(q",7") (22) 

+ (/3"+i/2)^VqD(q")p"+i/2 
+ ^(/3"+'/')^VqE(q",7")/3"+i/2 

/3"+i/2 = /3"-^(v^C/(q",7") + i(/3"+i/2)Tv^E(q",7")/3"+i/2) (23) 

*For the fixed lapse system H15|l we ignore the dependence on the discrete densitized lapse 
in the notation. 
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q" + _E((S(q"+l) + S(q"))p"+l/2 (24) 
+ (D(q"+l)^ + D(q")^)/3"+l/2) 

7" + ^((E(q"+i,7"+i) +E(q",7"))/3"+i/2 (25) 

+ (D(q"+l)+D(q"))p"+l/2) 
p«+l/2 _ ^(^l(pn+l/2)Ty^g(q„+l)p„+l/2 ^ ^^^^q„+i^^„+i^ 

+ (/3"+l/2)TVqD(q"+l)p"+l/2 (26) 
+ 1 (^«+i/2)T^^E(q"+l, 7"+l)/3"+l/2) 

/3"+'/'-^(^(V^C/(q"+\7"+') (27) 
+ i(/3»+i/2)T^_^E(q"+i,7"+i)/3"+i/2y 

In this system the computationally most expensive terms are subsumed in 
the expressions VqJ7 and "V^U . There is only one evaluation per step of those 
terms. 

The system (P^ - ([?71) ean be decomposed into three substeps, the first two, 
((I22l),(l23|)) and {^,^) being implicit in (p"+i/2, ^"+i/2) (qn+i .^n+i-,^ 
respectively. They are solved by fixed-point iteration, which is local at every 
grid point. The last substep f f^ . l^T)) ) is explicit. 

With that method the constraints are not explicitly taken into account. Yet, 
since in the continuous equations of motion the constraints are satisfied for all 
times as long as they are satisfied for the initial data, one may expect that the 
constraint violations in the discrete case remain small. 

3.2 The RATTLE method 

The RATTLE method ([1], [HI Section VII. 1], [H Chapter 7]) is an exten- 
sion of the Stormer-Verlet method to holonomically constrained systems. It is 
symplectic and time-reversible, of second order accuracy, and enforces both the 
holonomic and the derived hidden constraints in the numerical solution. When 
applied to (fT8|) .(ll9 |) . a step of the RATTLE method consists of the following 
equations, which form a nonlinear system for q"+^, 7"^^, p"^^ and (3"^^. 

1. First half-step for the momentum variables: 

P" - ^(i(p"+'/')'^VqS(q")p"+l/2 + VqC/(q",7") 

+ (/3"+i/2)7'VqE(q",7")/3"+i/2 (28a) 
+ (/3"+i/2)^VqD(q")p"+i/2 + Vqg(q")A"'+) 

/3" - ^((/3"+i/2)T^_^E(q",7")/3"+i/2 

+ V^[/(q",7")+At"^+) (28b) 



^n+l 
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n+1 _ 



,n+l 
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n+1 



n+1/2 



/3 



n+1/2 _ 
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2. Full step for the position variables: 

qn+i ^ q" + ^((S(q"+i) + S(q"))p"+i/2 (29a) 
+ (D(q"+l)^ + D(q"f)/3"+l/2) 

7"+' - 7" + ^((E(q"+\7"+i) + E(q",7"))/3"+i/2 (29b) 
+ (D(q"+i) + D(q"))p"+i/2) (7" ^ 0) 

3. Position constraints: 

g(q"+i) = (30a) 
7"+i = (30b) 

4. Second half-step for the momentum variables: 

-^(i(p"+'/'rVqS(q"+l)p"+V2 + VqC/(q"+\7"+^) 
+ (/3"+i/2)^VqE(q"+\7"+i);3"+i/2 (3^j^) 

+ (/3"+l/2)7^VqD(q"+l)p"+l/2 + Vqg(q"+l)A"+l^-) 

+ V^C/(q"+\7"+i) + M"+i'-) (31b) 

5. Momentum constraints: 

Vqg(q"+')^S(q"+i)p"+i + A(q"+i)/3"+i = (32a) 

D(q"+i)p"+i = 0. (32b) 

In the last equation (I32bp we used the fact that for our Hamiltonian the matrix 
E(q, 7) vanishes when 7 = 0. 

Equations (l28l)-(l3Ql) determine q"+\ and the system (l3T|) , ([321) determines 
pTC+i^ The equations can be solved by an iterative procedure that requires only 
the solution of linear systems with the matrices A(q") and A(q"+^) and their 
transposes. This procedure is described in 28], and we do not discuss it here. 

4 1+1 dimensional test cases 

In this article we are interested in the numerical properties of the discretized sys- 
tems ([5]) and (|15p . Our numerical experiments are based on a simplified model 
that arises in certain highly symmetric space-times. Essentially we assume that 
the 3-metric is diagonal, depends only on the coordinate and has just two 
independent components, — Qiii- The details are described in |28j. 

We distinguish between two essentially 1+1 dimensional classes of solutions 
of Einstein's equations, namely the spherically symmetric space-times where 



10 



C '■— sin^ and a second class where C = 1- The latter includes a perturbed 
Minkowski geometry. 

It is then natural to define 



i {h22 + C^h^s) 



7r22 + (tt 



33 



(33) 



and it turns out that fr is indeed the canonical momentum corresponding to h. 

In the spherically symmetric case we now consider the equatorial hyper- 
surface, i.e., x^ — tt/2 and ^ = 1. Using the new variables h and tt one can 
derive 1+1 dimensional counterparts for the Hamiltonians Tis, Ti-g and Ti.f3. The 
"1+1-dimensional Super-Hamiltonian" becomes 

-TT^^TT^^hiihii — TT^^irhiih 



nl ^ / dx 



a -dhdh - 2hd^h + hdhd\og{hii) + 2^hiih 



+ 27r"/iiia/3 + 7r"/?5/iii + n/Sdh 



and the 1+1-dimensional gauge Hamiltonians are 

1 



■Hi 



Hp 



dx 



= dx 



2a^h-fdh + ah^jda - -a-^h^Y - PldP 



(3ada — aad(3 + 2^'^'^^ 



- ~ 2 ~ 2 ~ 

a'^hjdh+ -a^h^-fd\oghii + -ah^^da 



(34) 



(35) 



(36) 
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Here we have ^ = 1 in the spherically symmetric case and ^ = if C = 1 • 

These Hamiltonians have a similar structure as their 3+1-dimensional coun- 
terparts (HI) and ((5|). p^ . their discretization will therefore be of the form (fT7)) . 



4.1 Space Discretization in the 1+1 dimensional setting 

The spatial discretization of the Hamiltonians ([34|) . ([35|) and ([36|) is done simi- 
larly as in |28| . We use piecewise constants functions to approximate /in, ft,, a, 
7, TT^^, TT, a and /3, where f3 and 7 are discretized on a staggered grid. Yet, for 
the discrete spatial derivatives we use a modified approach. 

We distinguish between four types of terms. The first three types contain 
only one (first or second) derivative, they are discretizations of terms of the form 
f^PdP or f^pd'^p, with functions /^, p and p. The fourth type contains 
a product of first derivatives, it is the discrete counterpart of f^dpdp. 

1. In the first type of terms the discrete counterpart of the functions /' are 
defined on the same grid. An example is the term 2ahd^h in ([M]) . For 
these terms we approximate the derivative with centered second order 
finite differences as in pS] : 

dp{x,) ^ (DoPh = (/iVi - /ti)/2Ax, 
d'Pix.) - (D+D^p), = {fl, 2p + p_,)/Ax\ (37) 
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2. In the second type of terms the discrete function /'^ is defined on the 
staggered grid and on the non staggered grid. An example is the 
term 27r^^ft,ii9/3 in ((M)) . In this term we discretize the derivative using 
one-sided finite differences 

{ffdf){x,) ^ flff{D^f% = flffUf - fl,)/Ax. (38) 

Since /"^ is defined on the staggered grid this formula is second order 
accurate at the grid points of the non staggered grid. 

3. For the third type of terms the discrete functions and are defined 
on the non staggered grid and on the staggered grid. An example is 
the term Tr^^f3dhii in p4p . For this term we average the function and 
use one-sided finite differences for the derivative of f^: 

iffdf){x.) - {fl+fl,)ff{D+f%/2 = 

= {fl + fl+i) fhf!+i - /f )/2Ax. (39) 

Again, since is defined on the staggered grid this formula is second 
order accurate at the grid points of the staggered grid. 

4. Finally, the fourth type of terms deals with discretizations of f^df^df^, 
where the functions /* are defined on the same grid. An example is the 
term adhdh in (|34p . There we use a combination of one-sided finite dif- 
ferences 

{fdfdf){x,) ^ fl{{D^f UD^f% + {D+f UD+f%)l2. (40) 
Again this formula is second order accurate. 

There are two reasons to change the discretization of spatial derivatives. At 
first, in [28j we observed instabilities with periodic boundary conditions and an 
even number of grid points. 

We can explain these instabilities as follows. When we discretize using cen- 
tered finite differences, i.e. f^dpdp{xi) — *■ fl{DoP)i{Dof'^)i, then the deriva- 
tive of the resulting Hamiltonian with respect to ff contains a term //(-Dq/^)* 
(see also e.g. [13] for the definition of Do, _D+ and _D_). This is a centered finite 
difference formula with an extended stencil {Dq discretization). If we linearize 
the corresponding equations and check the well-posedness of the result using 
the techniques proposed in [13] we find instabilities when the number of grid 
points is even. This is not the case when we discretize the spatial derivatives 
using pO|) . 

Moreover, when the discretizations ([37|) and (|40|l are used then we can (for 
periodic boundary conditions) perform partial integrations also in the discrete 
Hamiltonian, because 

i Y,f^ {{D+gUD+h), + {D^gUD^h),) = (41) 
= ^Y.^^a^{D+D^h), -\Y,g. {{D+JUD+h), + [D^fUD^h),) . 

i i 
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Boundary conditions and discrete gauge conditions. Concerning the 
treatment of boundary conditions we proceed analogously as in [28 . We apply 
periodic boundary conditions for the perturbed Minkowski problem and Dirich- 
let boundary conditions in the Schwarzschild space-time. The discrete boundary 
conditions are treated using ghost zones. 

As in [5H] the gauge conditions for the RATTLE scheme are derived from 
continuous gauge conditions of the form 

=0. (42) 

When we are interested in the generalized harmonic system (corresponding 
to Tig + Hg) then we choose 

g{hii,h, a) = d [huh-^ exp{2MC{R - 1))) , (43) 

where ^ = 1 in the spherically symmetric case and ^ = when C = 1- When we 
perform calculations with the fixed lapse system (corresponding to Hg + Hp) 
then we set 

gihii,h,a)=d[ahiih'^y (44) 

The reason for this choice is clearly that the analytic solutions for which we test 
the schemes (see section satisfy those gauge conditions. 

For the generalized harmonic system in the spherically symmetric case (^ = 
1) we see that the gauge condition depends on the mass parameter M. Hence, 
if we did not know the analytical solutions then it was not possible to apply 
this gauge condition. However, in principle one can also apply other gauge 
conditions, e.g. Dirac gauge [HEI, in the RATTLE method, but the results 
are then not compareable to the results of the free evolution scheme, because 
one needs to use different initial data (those that satisfy the gauge condition). 

4.2 Test scenarios 

In sections [5H3 we perform numerical experiments with the Stormer-Verlet and 
the RATTLE method. There we check how well the numerical schemes repro- 
duce a perturbed Minkowski space-time and the Schwarzschild space-time in 
two different coordinate systems. The analytical solutions are the following. 

4.2.1 A perturbed Minkowski metric 

The Minkowski metric describes a flat space-time, the analytical solution is 
(with X — x^, y — x'^, z — x^) 

ds^ ^ -d£- + dx^ + dy^ + dz"^ . (45) 

It is easy to check that for t — const, slicing this solution really is in the class of 
solutions with C = 1 and thus the numerical schemes we described are applicable. 
Here we perturb the Minkowski initial data such that 

/ill = 1 -f' £11, /i=l+e, a^l+Ea, 7 = £7 (46) 
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and 



n^-^S'', n = S, cr = S„, P = 6/3. (47) 

In the numerical examples we will choose the perturbations e, (5 to be Gaussian 
functions of width 1/20 and height 10^^. To avoid problems with boimdaries 
we apply periodic boundary conditions. 

4.2.2 Schwarzschild space-time 

In the numerical experiments that deal with the Schwarzschild space-time we 
consider two different coordinate systems. They are chosen such that we get 
solutions of the equations of motion for the Hamiltonians Hg + Hg and Hg + 
respectively. We denote these coordinate systems "radially harmonic" and 
"fixed lapse coordinates" respectively. 

Radially harmonic coordinates. From the Schwarzschild space-time in stan- 
dard coordinates (i, r, 6, (f) we come to the radially harmonic coordinates via the 
coordinate transformation 

^(^)^1+2^1'^K7^)' -(^^- l-exp(2M(l-i^)) - 

That is, R G (1, 00) and i? = 00 is the horizon, whereas i? = 1 is spatial infinity. 
The 4-metric in the new coordinate system is 

where dfi^ is the volume element of the sphere and r = r{R). From the 4-metric 
we read off the following expressions 

/ill = r^e^^'^fi"^), h^r\ a ^ r-\ 

tt" = 0, ^i- = 0, (3 = 0. (50) 

This solution satisfies the canonical Hamiltonian equations of motion that cor- 
respond to + TOg and the gauge condition for ^ = 1. 

Fixed lapse coordinates. The fixed lapse coordinates for the Schwarzschild 
space-time are obtained from the standard Schwarzschild coordinates (t, 0, </>) 
via the coordinate transformation 

i?(0 = (^)", r(i?) = 2Mi?Vii. (51) 

The 4-metric in the new coordinate system is 
ds^ = (l - i?-/") + ,,,^20/rf(T-i^-Vn) -^^^ + 4M^i^Vn,^2 
and we get 

4M2 11 _ 

" " 121i?20/ii (1 _ ' ^ -0, 

/i = 4M2i?2/i\ ^^0, (53) 

a = ll(2M)-3i?8/ii ( 1 _ R-^nA , p = o. 
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Figure 1: Perturbed Minkowski space-time: The maximum norm of Hamil- 
ton and momentum constraints, C and Ci respectively. Left: Results of the 
Stormer-Verlet method. Right: Results of the RATTLE method. 



This solution satisfies the canonical Hamiltonian equations of motion that cor- 
respond to Hg + Tip and the gauge condition for ^ = 1. 

5 A perturbed Minkowski problem: comparison 
of the four schemes 

In this section we compare the results of the four numerical schemes that we 
obtain when we combine the 1+1 dimensional generalized harmonic system (cor- 
responding to T-Cg + Hg) and the fixed lapse system (corresponding to T-Cg + Hp) 
with the Stormer-Verlet and the RATTLE integrator. 

Simulation data. We apply the schemes to the perturbed Minkowski problem 
described in section 14.2.11 The perturbations are Gaussian functions of width 
1/20 and height e = 10"^. 

When we apply the RATTLE scheme then we set = 0. The justification 
for this choice is that in the constrained scheme the equation 7 = is enforced 
for each time step and for the initial data it is always easy to satisfy it. 

Moreover, we set £0, = for the fixed lapse system, because there a is an 
external field and we assume that the correct function is given in advance. 

In the simulations we use 50 grid points and the size of the time step is a 
quarter of the spatial grid spacing. At — Ax /A. 

Simulation results. We see in figured] that the results for both systems of 
equations are similar when the free integration scheme using Stormer-Verlet 
is applied. Moreover the RATTLE scheme applied to the 1-1-1 dimensional 
generalized harmonic system also provides similar results. 

In those cases the evolution is stable at least until t = 1000 and we do not 
see growing constraints or errors. In particular the Hamilton and momentum 
constraint functions stay at almost the same size of about 2 • 10-^^-2 • lO^^. 
However, these functions do not vanish. 

The RATTLE scheme applied to the 1-1-1 dimensional fixed lapse system 
behaves differently. By construction it is clear that the momentum constraints 
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become small. They are the hidden constraints of 7 = and in our imple- 
mentation of the RATTLE scheme they are enforced with an accuracy of at 
least 10"^". However, at some point the evolution cannot be continued in this 
scheme, because the algorithm to solve the nonlinear system does not converge 
within lO'' iteration steps. The problem there is to satisfy 7 = 0, i.e. (I30bp . 

Discussion of results. The results that we get for this example indicate that 
strong hyperbolicity of the equations of motion is still essential when we use 
symplectic integrators for Hamiltonian systems. 

In [3^ we discussed tests of the same problem. There we used discretized 
ADM equation in the evolution and got rapidly growing high frequency errors. 
This is not the case here, and at least for three schemes all other desireable 
properties that the symplectic integration of the ADM equations showed (in 
particular the conservation of the harmonic energies) are still valid. This behav- 
ior was expected in the beginning, because the initial value problem of strongly 
hyperbolic systems is well-posed whereas it may be ill-posed for systems that 
are only weakly hyperbolic (the ADM equations are weakly hyperbolic but not 
strongly hyperbolic) [191 [El ■ 

Now, out of the four schemes that we investigated there are three that pro- 
vide the expected results. From the fourth scheme, the constrained integration 
of the fixed lapse system, we see that to obtain favorable results it is not suffi- 
cient to start from a strongly hyperbolic system. 

From a computational point of view we find that the reason for the problems 
in that scheme is basically that the momentum constraints cannot be enforced 
with the required accuracy, because the algorithm to solve the nonlinear system 
([^ - ([50]) does not converge. 

The solution of this scheme itself points to the interpretation that singular- 
ities occur during the evolution. In particular we observe large gradients in tt 
and (3. It is evident that these problems are not purely numerical, since we get 
quite similar results when we take a smaller time step. At = Aa:/8 and also 
with higher resolutions, N = 100 and N = 200. 

It is not clear whether this is a problem of the coordinate system or of the 
underlying solution. The algorithm that we apply projects onto the (momen- 
tum) constrained hypersurface in each time step, and in particular in the first 
one. As the constraints are quite large for the considered initial data the re- 
sult of this projection can be that the problem is no longer Minkowski with 
a small perturbation, but a solution that develops a physical singularity after 
some time. However, when we start from initial data with smaller constraints, 
e.g. by taking the height of the Gaussian perturbations as e = 10~^, we can 
evolve to later times. 

Since the fourth scheme provides better results when we start from initial 
data with small constraints one may expect that constrained integration is fea- 
sible when the initial data are chosen appropriately. This is supported by the 
calculations presented in section [7| where we find that it can be even beneficial 
to use this scheme. 
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Figure 2: Schwarzschild space-time: The maximum norm of the Hamilton and 
momentum constraints, C and Ci respectively, in the generalized harmonic 
evolution system. Left: Results of the free evolution scheme. Right: Results of 
the constrained evolution scheme. 

6 Schwarzschild space-time: free vs. constrained 
generaUzed harmonic evolution 

Here we discuss the Stormer-Verlet and the RATTLE scheme for the 1-1-1 dimen- 
sional generalized harmonic system (corresponding to Hg+Hg). We apply these 
schemes to the Schwarzschild space-time in radially harmonic coordinates (see 
section 14.2. 2p . The mass parameter is chosen to be M = 1 and the boundaries 
are at i?; = 2 and Rr = 3. 



Simulation results. In figure [2] we see that for both, free and constrained 
evolution the Hamilton and momentum constraints are growing exponentially. 
The same behavior can be observed for the error of the functions themselves. 
At some point the evolution breaks down. 

This happens almost at the same time for both integration methods. But 
one can evolve to later times if one uses a finer grid. Moreover, with better 
resolution errors become quadratically smaller. 



Discussion of results. For the generalized harmonic system in the Schwarz- 
schild space-time the results of both schemes are qualitatively similar. This 
is in contrast to [H] where we observed that for the Schwarzschild space-time 
constrained symplectic evolution of ADM equations leads to better results than 
free evolution. 

Yet, the difference there was that in the constrained evolution we were able 
to enforce the momentum constraints. Here we can only enforce a combination 
of the momentum constraints and cr = 0. Hence, the boundary condition can 
lead to violations of the momentum constraints. 

Again the results here are better than with the ADM evolution presented 
in [5S]. There we found that the evolution breaks down at t « 20 for the 
constrained scheme and at t « 3 in the free evolution. Here we can evolve until 
<: « 40 with both schemes. Moreover, in contrast to the free evolution using the 
ADM equations, errors are smaller in a fine grid and we can evolve longer than 
in a coarse grid. 
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Figure 3: Schwarzschild space-time: The maximum norm of the Hamilton and 
momentum constraints, C and Ci respectively, in the fixed lapse evolution sys- 
tem. Left: Results of the free evolution scheme. Right: Results of the con- 
strained evolution scheme. 
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Figure 4: Schwarzschild space-time: The maximum of the error of /in, in the 
fixed lapse evolution system. Left: Results of the free evolution scheme. Right: 
Results of the constrained evolution scheme. 

7 Schwarzschild space-time: free vs. constrained 
evolution with fixed densitized lapse 

Finally we compare the Stormer-Verlet and the RATTLE scheme for the fixed 
lapse system (corresponding to Ti^ + 7i^). We apply these schemes to the 
Schwarzschild space-time in fixed lapse coordinates (see section [4.2.2p . Again 
the mass parameter is chosen to be M — 1 and the boundaries are at i?; = 2 
and Rr = 3. 



Simulation results. We see from figures [3] and S] that with the free evolution 
scheme we obtain in the beginning (until t k, 40) similar results as with the 
generalized harmonic system (see section[5|) . In particular constraints and errors 
are growing exponentially and become quadratically smaller when the resolution 
is increased. 

At later times errors are still growing exponentially, but it may happen that 
they are bigger when a better resolution is used and that we can evolve to later 
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times in a coarser grid. 

In the constrained scheme we find that at least until t = 200 the Hamilton 
constraint stays almost at the same size. For the errors of the functions them- 
selves one finds at most linear growth. The momentum constraints are of course 
very small by construction. 

Discussion of results. Here, as in [23 where we investigate the ADM system, 
we find better results when we apply the constrained RATTLE scheme. 

In the free evolution we observe exponentially growing errors and instabili- 
ties. Presumably the reason for these problems are the naive Dirichlet boundary 
conditions that lead to constraint violations. 

But in the RATTLE scheme there is no sign of growing constraints and 
moreover the errors of the functions themselves increase very slowly. Thus, for 
the systems and examples investigated so far constrained symplectic evolution of 
systems with boundaries leads to better results when the momentum constraints 
are enforced. It will be interesting to see whether this observation can be made 
in higher dimension, too. 

In 1-1-1 dimensions, as there are no gravitational waves possible, we only 
need to worry about the boundary conditions to respect the constraints. In the 
constrained RATTLE scheme this is only partially ensured, it is in particular 
surprising that there are no increasing violations of the Hamilton constraints. 

8 Conclusion 

This article deals with symplectic numerical integration of strongly hyperbolic 
Hamiltonian formulations of general relativity. We discussed two appropriate 
Hamiltonians that are constructed by introducing hyperbolic drivers for the 
shift and in one system also for the densitized lapse. These are the fixed lapse 
system and a generalized harmonic system derived by Brown [11| . respectively. 

For both systems we performed numerical experiments using a free (Stormer- 
Verlet) and a constrained symplectic integrator (RATTLE). The free evolution 
schemes lead to similar results for both systems. In contrast to the ADM evo- 
lution [28] we do not see rapidly growing high frequency errors. Also the good 
propagation properties of nearly conserved quantities that were found for some 
situations in the ADM evolution are recovered here. We thus conclude that it 
is indeed better to apply free symplectic integrators with strongly hyperbolic 
Hamiltonian systems. 

With the constrained schemes the situation is different. When we apply the 
RATTLE method with the generalized harmonic system then we always got 
results that are very similar to the free evolution schemes. Since constrained 
integration requires much more computational resources than free evolution we 
thus conclude that constrained integration is not practicable with the general- 
ized harmonic system. 

Concerning constrained integration with the fixed lapse system we have 
shown one example where we observe instabilities in the evolution, but also 
one example where this scheme is much more stable than the others. 

In the former case we argued that the problems occur due to constraint 
violations in the initial data. At the moment we cannot be certain that this is 
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really the case, but the results can be interpreted in that way. We move this 
question to future investigations. 

Yet, if the constraint violations in the initial data are indeed the reason for 
the problems then this RATTLE scheme is an interesting option to deal with 
the constraints for situations with appropriate initial data. This statement is 
supported by the results of that scheme for the fixed lapse and the ADM system 
(see [35]) with spherically symmetric solutions. There we observe stabilizing ef- 
fects in the constrained evolution schemes. Since we do not find these effects 
when we apply the RATTLE method with the generalized harmonic system we 
suppose that our constrained evolution scheme has advantages when the mo- 
mentum constraints are enforced in the equations of motion and the initial data 
satisfy the constraints!! It will be interesting to see whether this conjecture can 
be supported by analytical studies of the underlying elliptic-hyperbolic system 
of equations or by further numerical investigations. 

Another question to address in the future is the gauge choice. The systems 
that we consider here are quite restrictive for that matter. In particular in 
the fixed lapse system the requirement of strong hyperbolicity already fixed the 
Hamiltonian completely. This forced us to use special coordinates in the test 
examples, a factor that is clearly not satisfactory. One question is thus whether 
appropriate Hamiltonians exist that allow e.g. general Bona-Masso slicings. 

We also want to start to investigate examples in higher dimensions. We 
discussed that the free symplectic integrators are implicit, but that this im- 
plicitness is very mild such that symplectic integration does not need more 
computational resources than standard explicit methods. However, symplectic 
integrators should be applied to Hamiltonian systems. We expect that the right 
hand sides of discrete systems with that property will be more complicated than 
non Hamiltonian ones. But the size of the difference is not yet known. 
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A Derivation of the fixed lapse Hamiltonian 

In this article we dealt with two strongly hyperbolic Hamiltonian formulations 
of general relativity, namely the generalized harmonic system (jSl) and the fixed 
lapse system (ITSI) . The former is similar to the generalized harmonic system 
derived in Here we discuss the derivation of the latter. 

As explained in section lT^ for this system the densitized lapse is an external 
field and the shift becomes a momentum variable. The canonical symplectic 
two- form for the phase-space of (/ly , 7^; tt*-' , is then dhij A dTr*-' -|- d% A d/?*. 

Furthermore we assume the following form of the Hamiltonian 



^The momentum constraints are enforced in the RATTLE scheme for the ADM equations 
and the fixed lapse system, but not for the generalized harmonic system. 




(54) 
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with a function that needs to be chosen appropriately. 

To obtain this function we naturaUy start from the ansatz that Brown pro- 
posed in [llj . Yet, here we are interested in the principal part of the resultinj 
equations of motion only and omit all terms that do not contribute there 
Hence, we start with 



-f3-'djf3 — 02'(.a ^ jk"-" ~ ^3 



(55) 



where C2, C3, C5 and Cg are constant parameters. 

Having this ansatz for Cf^ we can derive the equations of motion that corre- 
spond to the Hamiltonian (|54p and their principal part. Wc obtain 
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where the symbol = is used to denote equality up to lower order terms and for 
the non vanishing coefficients in this principal part we get 



\ klm 

Aim 

ijk 



A 



iklm 



^.^UUikf^lm 



lie. 



A im 

Cijki 

^ik 

jjijklmn 



SY'hjk + SJ'hik, 

Sir, 

= 2ahikhji - ahijhki, 



C3 + C5) a^hU'^'h^^ 



'-'kl 



J^kn^tmyl _ hlnf^tkf^jm _^ hmnf^tkf^jl ^ 2/l'^'/l™/l^" 

{C2~C3 + C5)h'^h 



k7n\ 



Gm 
ikl 

^km 



2 

-2hik5l 



(56) 



The remaining coefficients C", C, D' , D and D vanish. 

Now, according to J19' this second order in space evolution system is strongly 
hyperbolic if the second order principal symbol 
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®In the terms that do not contribute to the principal part are needed in order to 
ensure that becomes a covariant vector under spatial diffoomorphisms and a scalar density 
of weight +2 under time reparametrizations. 
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has real eigenvalues for all unit vectors and is diagonalizable in a regular 
way. Therefore the aim is now to choose the parameters C2 , . . . , Ce such that 
P becomes diagonalizable. 

The difference in comparison to the analysis of other formulations, like e.g. 
NOR [ZF , is that we do not get strong hyperbolicity for a whole range of pa- 
rameters. Instead we can derive nonlinear relations that the parameters need to 
satisfy. That is, only on a lower dimensional subset of the parameter space we 
obtain appropriate formulations. It is then very helpful if one has an algorithm 
to derive these relations. We apply the following one. 

A.l Diagonalizability of a parameterized matrix 

Starting from a matrix V that depends on parameters Ci we want to choose the 
parameters such that V is diagonalizable. The idea to find appropriate param- 
eters is to adapt the algorithm for the calculation of the Jordan decomposition 
of a matrix. 

The first step is to calculate the eigenvalues of V, i.e. to find the roots of 
det{V — A I). When V has a general form already this step might fail, because 
the analytical expressions for the eigenvalues are needed and for quasilinear 
systems the components of V depend on the fields. Moreover, the degree of 
the polynomial det(7' — A I) can become big and one obtains long expressions 
soon. However, for our problem we indeed get analytical expressions for the 
eigenvalues. 

Now, the multiplicity of the eigenvalues of V can be one or bigger than one. 
For the former eigenvalues nothing needs to be done, but for the latter it may 
happen that there are less eigenvectors than the multiplicity of the eigenvalue. 
In this case we proceed as follows. 

Let Ai be an eigenvalue of multiplicity two or more. We calculate the kernels 
Ki = ker(V - Ai I) and K2 = ker{V - Ai I)^. If Ki = K2 then the geometric 
multiplicity of Ai is the same as its algebraic multiplicity and nothing needs 
to be done. But if there are vectors in K2 that are not in Ki then the set of 
eigenvectors of Ai is not complete. That is, we need to restrict the possible 
choices of the parameters Ci. 

To derive the corresponding equations we pick some vector x E K2 \ Ki and 
calculate y = {V — Xi l)x. The vector y depends on the parameters and we 
probably can choose them such that ?/ = 0. This leads to the desired relations. 
Of course these relations must not depend on the fields. It may happen that we 
cannot choose the parameters appropriately. In this case we don't get a strongly 
hyperbolic formulation. 

Finally we use the derived relations to reduce the number of parameters in 
V and eventually repeat the procedure if there are still eigenvalues for which 
the geometric multiplicity is smaller than the algebraic one. 

A. 2 The fixed lapse system 

Now, to analyse the diagonalizability of P (defined in ((57)) ) with the algorithm 
described in the previous section lA.ll it is helpful to introduce an orthonormal 
basis {n*, w*, w*}. When we expand the tensor indices in hij, tt*^ , /3* and 7^ into 
this basis, denoting e.g. h^y = hijV^v^ , then P decomposes into three blocks. 
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The first block, Pi, corresponds to the subsystem of the components hy^u 
and tt"'". Its eigenvalues are /3" ± a^/h and it has a complete set of eigenvectors 
for each choice of the parameters C2, ■ ■ ■ ,Cq. 

The second block, P2, comes from the subsystem of the components hnm hyy, 
hww, tt"", tt"", tt™™, /?" and 7„. Its eigenvalues are /3" ± aVh with multiplicity 
one and /3" ± a^h{l + C2 + C3 — C5) with multiplicity two. For a general 
choice of the parameters it turns out that P2 is not diagonalizable. 

We first consider the eigenvalue Ai = If we apply the algorithm described 
in section lA.ll then we find that diagonalizability of P2 implies 

Ce^-l (9C2 + IOC2C3 + Cl - IOC2C5 - 2C3C5 + Cl) . (58) 


Then, replacing Cg in P2 using ([58| and applying the same steps as before for 
A2 = /?" - a^h{l + C2 + C3- C5) we get 

C2 = ^(4-C3 + C5). (59) 

Here we assume that the parameters are chosen such that there are no addi- 
tional degeneracies of eigenvalues. In particular we must choose C2 + C3 — C5 ^ 
— 1 and C2 + C3 — C5 7^ 0. Furthermore, the reality condition on the eigenval- 
ues implies C2 -I- C3 — C5 > —1. Now, choosing C2 and Cg in this way P2 is 
diagonalizable with real eigenvalues. 

What remains is to analyse the third block of P, namely P3. It comes from 
the subsystem of hnv, h„w, tt"", tt"*^, 7„ and ^wl\ When we use ([55)) 

and ([5^ to remove C2 and Cg from P3 then we obtain the eigenvalues /3" ± 
l/3a^h(4: — C3 + C5), both with multiplicity four. The algorithm of section 
lA.ll applied with any of these eigenvalues then leads to 

C3 = -^ + C5. (60) 

Hence, altogether we get 

C2 = ^, C3 = -y + C5, C6 = ^. (61) 

If these relations are satisfied then the eigenvalues of P are automatically real. 
Hence, we have a strongly hyperbolic system and for C5 = 2/7 we obtain the 
fixed lapse system ([T5|). 

Since P only depends on the difference C3 — C5 and not on the sum of 
these parameters it is clear that other strongly hyperbolic formulations that are 
derived from the ansatz (|55p have the same principal part. There might be an 
exception when the parameters are chosen such that some eigenvalues coincide 
(e.g. when C2 + C3 ~ C5 ~ 0). However, for our purposes it was sufficient to 
find a single strongly hyperbolic Hamiltonian system, and here we were able to 
derive one. 

The calculations described above were performed using Mathematica [30j 
and the open-source package xTensor for abstract tensor calculations developed 
by J. M. Martm-Garcfa [H]. 

''The first and third block do not appear for the simphfied Hamiltonian T-Og + because 
the corresponding system is just 1+1 dimensional. 
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